Reversible self-assembly of patchy particles into monodisperse icosahedral clusters 
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We systematically study the design of simple patchy sphere models that reversibly self-assemble 
into monodisperse icosahedral clusters. We find that the optimal patch width is a compromise 
between structural specificity (the patches must be narrow enough to energetically select the desired 
clusters) and kinetic accessibility (they must be sufficiently wide to avoid kinetic traps). Similarly, 
for good yields the temperature must be low enough for the clusters to be thermodynamically 
stable, but the clusters must also have enough thermal energy to allow incorrectly formed bonds to 
be broken. Ordered clusters can form through a number of different dynamic pathways, including 
direct nucleation and indirect pathways involving large disordered intermediates. The latter pathway 
is related to a reentrant liquid-to-gas transition that occurs for intermediate patch widths upon 
lowering the temperature. We also find that the assembly process is robust to inaccurate patch 
placement up to a certain threshold, and that it is possible to replace the five discrete patches with 
a single ring patch with no significant loss in yield. 

PACS numbers: 81.16.Dn,47.57.-s 



I. INTRODUCTION 

The remarkable ability of biological matter to robustly 
self- assemble into well defined composite objects excites 
the imagination, suggesting that these processes could 
perhaps be emulated through the judicious design of syn- 
thetic building blocks. 1 Viruses provide a particularly in- 
spiring example. As first shown 40 years ago for the cow- 
pea chlorotic mottle virus^ and later for a wide variety 
of other species^^ 5 - empty spherical virus shells (cap- 
sids) can be made to reversibly self-assemble from indi- 
vidual protein subunits (capsomers) in vitro, simply by 
changing solution conditions such as the pH. Although 
this process resembles micellar self-assembly, there is an 
important difference: in contrast to the polydisperse dis- 
tributions that characterise micelles, the complete virus 
capsids are monodisperse. 

Of those viruses with a roughly spherical structure, the 
vast majority are found to have icosahedral symmetry. 
Many of the structures observed were first predicted on 
symmetry grounds by Caspar and Klug. 6 More recently, 
Zandi et alii and Glotzer et al4 found that the most 
stable sizes of simple model systems correlated with many 
of the structures found in virus capsids. 

By contrast, the mechanisms by which capsid 
self-assembly takes place are incompletely under- 
stood, though there has recently been consid- 
erable theoretical 9 -^^^^^^ and 
experimenta l 21 i 22 ' 23 i 24 ' 25 i 26 i 27 progress. Some com- 
mon conclusions emerge from this body of work. Firstly, 



one of the major potential obstacles to successful 
assembly results from a class of kinetic traps, in which 
a large number of partial capsids are formed but are 
then unable to proceed to completion because of a lack 
of remaining free capsomers. This situation can be 
avoided if the free energy barrier for nucleation is large, 
such that intermediates are formed only slowly and free 
capsomers remain plentifuL 16 ' 17 i 18 ' 19 i 20 ' 21 i 22 ' 23 Secondly, 
it has been shown that only weak capsomer-capsomer 
interactions are required for capsids to be stable, be- 
cause of the large number of interactions present in a 
complete capsidi 19 ' 20 i 21 ' 22 i 23 ' 24 It is worth noting that 
weak interactions will tend to result in a large nucleation 
barrier, and also allow easier breakdown of incorrectly 
formed intermediates. These factors together imply 
that in general weak interactions will favour successful 
assembly, so long as they remain strong enough for the 
complete capsid to be stable and for nucleation to occur 
in a reasonable time scale. 

Understanding how self assembly with virus-like con- 
trol and fidelity could be achieved with synthetically 
produced sub-units is an important goal for nanofab- 
rication. New experimental techniques to create self- 
assembling systems are being rapidly developed* 1 - For ex- 
ample, one particularly impressive recent example is the 
reversible formation of tetrahedra from four appropri- 
ately designed DNA strands 3? For colloidal and nanopar- 
ticle building blocks, one of the most crucial requirements 
to allow increased control over assembly is the synthe- 
sis of anisotropic "patchy" particles. There has recently 
been a great deal of work in this area, with notable sue- 



2 



cesses for both colloid s 29 ' 30 ' 31 ' 32 ' 33 and polymer-coated 
nanoparticles , 34 i 35 i2£ 

Both these experimental advances and the impor- 
tance of understanding capsid formation have stimulated 
a number of simulation studies of the self-assembly of 
monodisperse objects. Rapaport^ 4 - produced the first 
simulation study of capsid assembly, and used a model 
consisting of trapezoidal particles with attractive sites 
which assemble into a virus-like structure. However, Ra- 
paport's model uses several unphysical rules to assist in 
correct assembly, such as irreversible bonding if subunits 
come together in the correct conformation. 

Subsequently, there have been a number of simulations 
that have achieved reversible self-assembly 15 ' 18 i 37 ' 38 i 39 
For example, Hagan and Chandler were able to assem- 
ble large 60-particle virus models by using as their ba- 
sic "capsomer" units spherical particles with directional 
anisotropic interactions that are chosen to be comple- 
mentary in order to help guide the particles into the right 
local relative orientations This feature may naturally 
be realised in protein-protein interactions but may be 
more difficult to implement in synthetic systems. Van 
Workum and Douglas also found that virus-like assem- 
blies (not necessarily monodisperse) formed from parti- 
cles made up of three dipolar Stockmayer particles con- 
nected together to form triangles. 37 Zhang and Glotzer 
created anisotropic model particles by rigidly connecting 
spheres with different attractive potentials, and showed 
that these can form a rich variety of small cluster shapes 
and extended structures! 38 ' 39 Most recently Nguyen et al. 
performed extensive studies using a geometric model in 
which appropriately shaped particles with specific inter- 
actions reversibly assembled into icosahedral capsids^ 
They produced a phase diagram for the model, and 
also observed incorrectly formed "monster" particles very 
similar to those seen with real viruses. 

In this paper we describe a set of minimal models, sin- 
gle spheres with anisotropic "patchy" interactions, that 
exhibit reversible self-assembly into monodisperse clus- 
ter phases. Choosing such simple systems facilitates the 
systematic exploration of parameter space to uncover the 
optimal design rules for self-assembly, and helps untan- 
gle the roles of thermodynamic and kinetic factors. This 
model may offer insight into how synthetic nanoparticles 
and colloids could be designed to self-assemble, and while 
the model is too simplistic to accurately describe the in- 
teractions between the capsid proteins of viruses, we hope 
that the principles identified will also be of relevance to 
biological examples of monodisperse self-assembly. 
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but the attraction is modulated by an orientationally de- 
pendent term, V ang . Thus, the complete potential is 



Vhj(rij) r < <j l ,i 

VhJ (r-y)14ng (f ij ,fli,flj) r > cr LJ , 

where is the orientation of particle i. V ang has the 
form: 



^ang (j^ij , , ) — Grij {^ij , ) Gji (^ji , ) (3) 

G tf (fy.flO = expf ] , (4) 



where a is the standard deviation of the Gaussian, 9kij 
is the angle between patch vector k on particle i and 
the interparticle vector ry , and & m i n is the patch that 
minimizes the magnitude of this angle. Hence, only the 
patches on each particle that are closest to the interpar- 
ticle axis interact with each other. 

When using this potential in our Monte Carlo simula- 
tions, we truncate and shift the potential using a cutoff 
distance of 3<7lj in order to improve the computational 
efficiency. To maintain the isotropic nature of the repul- 
sion, we also adjust the distance at which we start to 
include the angular modulation so that it corresponds to 
where the cut-and-shifted Lennard- Jones potential passes 
through zero. 

The target structure we choose to try to assemble is the 
12-particle hollow icosahedron (Fig.QJb)). This structure 
is sufficiently complex to allow interesting behaviour to 
be observed, but forms easily enough for large numbers 
of simulations of acceptable length to be carried out. It 
is also analogous to the structure of "T = 1" viruses. 

In order to design the patchy particles so that they 
self-assemble into hollow icosahedra, the natural choice 
for the patch positions is to place them such that they 
point directly at the neighbouring particles in the target 
structure. The resulting particle design, along with an as- 
sembly of twelve such particles into the target structure, 
is shown in Fig.QJb). Having chosen the patch positions, 
the only parameter in the potential that remains to be 
optimized is then the patch width a. 



B. Monte Carlo simulation 



II. METHODS 
A. Model 

Our model consists of spherical particles with a number 
of patches whose geometry is specified by a set of patch 
vectors. The repulsion is based on the isotropic Lennard- 



We simulate the system in the canonical ensemble us- 
ing a Metropolis Monte Carlo algorithm where the al- 
lowed move types are translations and rotations of indi- 
vidual particles. The translational moves are randomly 
chosen from a small cube centred on the selected parti- 
cle. Rotational moves make use of a quaternion repre- 
sentation of the particle's orientation, which is modified 
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by addition of a smaller random quaternion and then 
renormalized. 40 

This algorithm does not directly simulate physical tra- 
jectories, and is more often used to obtain thermody- 
namic averages of properties of a system. However, while 
the dynamics generated are not strictly rigorous, the use 
of small translational and rotational steps will mimic dif- 
fusive dynamics. 41 ' 42 ! 43 The advantage of this technique 
is that the derivatives of the potential are not needed, 
which allows the simulations to run very efficiently. 

In the analysis of the results it is useful to assign the 
particles to clusters. The condition used is that two par- 
ticles are considered bonded if their interaction energy 
is less than — 0.4 e. Two particles are then considered to 
be in the same cluster if they are joined by a continuous 
chain of bonds. 

As the above Monte Carlo simulations generally do 
not reach equilibrium, to determine some of the thermo- 
dynamic properties we also used the umbrella sampling 
technique. 44 ! 45 In umbrella sampling rather than gen- 
erating configurations with the Boltzmann distribution, 
configurations are instead sampled according to the dis- 
tribution exp [— (V + w(Q))], where w(Q) is a weight- 
ing function that is a function of an order parameter 
Q. Generally, the reason for introducing the weighting 
function is to artificially reduce the free energy barrier 
between different (meta)stable states of a system, hence 
facilitating interchange between these states and allowing 
the system to approach equilibrium. Canonical averages 
are simple to obtain from a simulation of such a non- 
Boltzmann (nB) ensemble using the expression 

(5> WT = <£exp[/MQ)]) nB , (5) 

where B is some generic property of the system. 

In our case, we wish to locate the temperature at which 
the system becomes most stable as a collection of icosahe- 
dra, and so we need to overcome the free energy barrier 
for the formation of this structure. As the icosahedra 
only very weakly interact with each other, it is sufficient, 
to a first approximation, just to consider a 12-particle 
system. For an order parameter we simply use the num- 
ber of particles in the largest cluster present. We wish to 
choose the weight function so that the system spends an 
approximately equal amount of time at each value of Q. 
To achieve this we use an iterative scheme in which w(Q) 
is updated at the end of each of a series of simulations. 



III. RESULTS 

A. Energetics 

A minimal requirement for self-assembly is that the 
system of monodisperse clusters is energetically most sta- 
ble, i.e. a system of n icosahedra must be lower in en- 
ergy than the lowest-energy single 12n-particle cluster. 
We check for what range of a this requirement holds, by 
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FIG. 1: (Colour Online) (a) The size dependence of the en- 
ergy of the global minimum for different values of a. Each 
line is labelled by the value of a (measured in radians), with 
particularly stable cluster sizes also labelled. The energy zero 
is a fit to the energies of the assemblies of complete icosa- 
hedra possible at JV=12, 24, 36 and 48 at the appropriate 
value of a. Hence, negative values indicate that it will never 
be energetically favourable for the system to form monodis- 
perse icosahedra. (b) A single particle designed to form a 
hollow icosahedron, and a complete icosahedral cluster. The 
yellow spheres represent the positions of the patch vectors, 
(c) The minimum-energy structures for 24 and 36 particles 
for g < 0.6. The radius of the particles is reduced compared 
to (b) to allow the structures to be seen more clearly. 



performing global optimisation to find the lowest-energy 
configurations of our system as a function of the number 
of particles. The global optimisation algorithm we used 
was the basin hopping algorithm, 4e which has proved to 
be particularly successful for clusters4i Fig.[T]shows that 
for wider patches the global minima correspond to single 
compact clusters with magic numbers as for the isotropic 
Lennard-Jones potential (iV=13, 19, 23, . . . )4£ However, 
for sufficiently narrow patches (a < 0.6) the most stable 
sizes occur at ^=12, 24, 36, . . . , and do indeed corre- 
spond to packings of 1, 2, 3, ... discrete icosahedra (Fig. 
[He)), as desired. (Note that as the icosahedra weakly in- 
teract, the discrete icosahedra are touching in the global 
minimum configuration) . The energy landscapes for such 
magic number clusters are likely to have a single-funnel 
topography . 49 ! 50 



B. Self assembly 

To further investigate the self-assembly behaviour of 
our system, we performed Monte Carlo (MC) simulations 
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FIG. 2: (Colour Online) Snapshots of a system of 72 particles 
assembling into six complete icosahedra at a = 0.45 and T = 
O.Utkg 1 after (a) 3000, (b) 8000, (c) 60 000 and (d) 250 000 
MC cycles. 



in the canonical ensemble. Fig. [2] illustrates a typical 
trajectory, and shows that, given the right conditions, 
the patchy particles are able to reversibly assemble into 
a monodisperse set of hollow icosahedral clusters. At 
the beginning of the simulation, disordered clusters form 
rapidly (Fig. [Ha)). These clusters then soon begin to 
order. For example, five-coordinate vertices become vis- 
ible, and many of the particles already occupy their cor- 
rect positions in partially formed capsids (Fig.^b)). The 
slowest process is then the diffusion of the few remaining 
monomers to the empty sites in the clusters (Fig. [2Jc) 
and (d)). 

The simplicity of our model allows us to map out in 
detail the assembly behaviour, thus enabling the opti- 
mal patch width and temperature to be identified. Fig. [3J 
summarizes results from 480 simulations at different com- 
binations of a and T, where each simulation was started 
from a disordered configuration generated at high tem- 
perature. Although this plot is for a particular density 
and after a certain number of MC moves, the generic 
features are not sensitive to these parameters (for the 
dependence on density see Section IIIIE[) . 

The key feature of Fig.[3]is that there is a limited range 
in the space of patch width and temperature where self- 
assembly occurs efficiently. The maximum yield of 88% 
occurs at a = 0.45 and T = 0.14 ekg 1 . The general pat- 
terns shown in Fig. [3] are relatively insensitive to simula- 
tion time, although in general yields increase with time 
For example, we have observed yields as high as 98% after 
800 000 MC cycles. 

To understand the assembly kinetics summarized in 
Fig. [3J we first examine the effect of temperature at the 
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FIG. 3: (Colour Online) (a) The number of icosahedra 
formed, (b) the mean cluster size (averaged over particles) 
and (c) the average number of bonds per particle (counting 
each bond twice) after 80 000 MC cycles as a function of the 
patch width a (measured in radians) and the temperature for 
1200 particles at a number density of 0.15o"lj. The white 
lines show the thermodynamic transition temperature T c i ust 
for the transition from icosahedra to a gas of monomers (see 
Section [1113. 



optimal a. There is a clear maximum in the number of 
icosahedra formed, with the number dropping to zero at 
high and low temperature. The upper limit is thermo- 
dynamic because at high temperatures the stable phase 
of the system is a gas of monomers. The lower limit is 
kinetic and arises because at low temperatures the sys- 
tem does not have enough thermal energy to escape from 
incorrect arrangements of the particles, and so large low- 
density kinetic aggregates grow. These two limits are 



clearly visible from the mean cluster size plotted in Fig. 
Mb)- 

A snapshot of a kinetic aggregate is shown in Fig. BJa). 
This structure is typical of those formed at higher a and 
low temperature. The relatively wide patches allow sig- 
nificant rearrangement of the clusters, and it can be seen 
that there is a large degree of local order, with most par- 
ticles forming strong interactions through most or all of 
their five patches. The diameters of the "strings" in the 
network are roughly equal to the diameter of an ideal 
icosahedron. However, the low temperature prevents 
complete rearrangement to form icosahedra. For lower 
a and low T the frustration of the system is much more 
severe, and most particles only manage to form two or 
three strong interactions (Fig. 03c)). The resulting net- 
work does not coalesce into dense clusters like those in 
Fig. BJa), and the "threads" of the network are only one 
particle in diameter. At very low a very few interactions 
are formed at all, and only small clusters are found. 

Having a temperature window where ordering can oc- 
cur such as that observed here is common. For example, 
in the protein folding community it has been argued that 
increasing the ratio Tf/T g (Tj is the 'folding' tempera- 
ture below which the native state becomes most stable 
and T g is the 'glass transition' temperature below which 
the dynamics becomes too slow for folding to occur) en- 
hances the ability of a protein to fold^i Similarly, here 
we find that the optimal value of a occurs roughly where 
there is the largest difference between T c i ust , the tempera- 
ture below which the clusters become thermodynamically 
most stable, and T g , below which ordering is kinetically 
hindered. 

The effect of reducing a from its optimal value is to 
narrow the temperature window over which icosahedra 
form. Firstly, T c iust decreases as the patches become 
narrower, because the vibrational entropy of the icosahe- 
dral clusters decreases. Secondly, the potential becomes 
"stickier" as the patches narrow, and so the 'glass transi- 
tion' temperature T g below which the particles becomes 
trapped in large non-equilibrium aggregates initially in- 
creases. Consequently, for a < 0.1 icosahedra are virtu- 
ally never seen. 

Fig. [Q shows that beyond a certain cr, the icosahedra 
are no longer the most stable particle configuration. This 
is mirrored in Fig. |3{b) by the large disordered clusters 
that form for larger a. At these values of a and den- 
sity the system phase separates into one or more liquid 
droplets coexisting with a vapour. A typical configura- 
tion from this region of the phase diagram is shown in 
Fig. Eflb) . Note that the liquid droplets are not partic- 
ularly spherical. The system has a low surface tension 
because most patches for the particles on the surface of 
the droplets are pointing inwards. As a result large shape 
fluctuations in the droplets are observed. 

The interaction between the clustering transition and 
the formation of a liquid phase gives rise to a particularly 
intriguing form for the phase diagram of this system. For 
example, as the temperature is lowered for a ~ 0.55, 



5 




FIG. 4: (Colour Online) Typical aggregate structures formed 
at non-optimal temperatures and patch widths, (a) String- 
like kinetic aggregates formed at T — 0.04efc~ 1 and a = 0.6. 
(b) Thermodynamic liquid-like aggregates formed at T = 
O.lSefc- 1 and a = 0.65. 



the system first condenses from a monomeric gas to a 
bulk liquid-like phase. However, at lower temperature it 
forms a gas phase again, but now made up of weakly 
interacting icosahedral clusters. This reentrant lower 
liquid-to-vapour transition is driven by the lower energy 
of the clusters, which overcomes the higher entropy of the 
liquid. Preliminary parallel tempering simulations have 
confirmed this unusual scenario. The probable implica- 
tion is that there is a closed-loop liquid-vapour coexis- 
tence line in the phase diagram for this system. This will 
be explored further in future work. 
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C. Thermodynamics 

In order to obtain estimates for T c iust, the temperature 
below which the system is most stable as a collection of 
icosahedra, we made use of the umbrella sampling scheme 
described in Section III Bl The resulting heat capacity 
plots obtained for a selection of patch widths are shown 
in Fig. |Sfa). 

The peaks in the heat capacity signal the transition 
from a gas of monomers to a gas of icosahedra. As an 
aside, we should note that this transition is not a phase 
transition, because it is fundamentally a transition asso- 
ciated with a finite cluster, e.g. the heat capacity peaks 
will retain the same finite width, no matter how large the 
system is. 

The positions of the peaks provide estimates for Tciust, 
which arc shown superimposed on Fig.[3Ja). It is interest- 
ing to note that relatively little supercooling is required 
for the icosahedra to begin to assemble. This is confirmed 
in Fig. [5tc) where we compare the equilibrium probabil- 
ity of a particle being in an icosahedron to that obtained 
in our MC simulations of the self assembly. In particular, 
when the icosahedra first become stable as the tempera- 
ture is decreased and their equilibrium probability is still 
low, the yield in our self assembly simulations mirrors (to 
within statistical error) the equilibrium results. 

The free energy profiles, A(Q), can be obtained 
straightforwardly from the umbrella sampling simula- 
tions, through the relation: 



A(Q) =A-kTlnp(Q), 



(6) 



and a series of profiles are plotted for the optimum patch 
width in Fig. [5£b). The change in the relative stability 
of the icosahedra compared to the vapour as the transi- 
tion is crossed is apparent, and for sufficiently low tem- 
perature the profile becomes barrierless. At T c iust, the 
free energy barrier is of the order of 2 kT. This rela- 
tively low value suggests that overcoming the nucleation 
barrier is not rate-limiting, and helps to explain why the 
yield of icosahedra initially follows the equilibrium line as 
the temperature decreases from above T c i ust (Fig. EJc)). 
However, as the equilibrium population of icosahedra in- 
creases, the yield from our self-assembly simulations be- 
gins to deviate from the equilibrium line, initially because 
the time scale for diffusion of free particles to the correct 
sites on incomplete clusters begins to limit the yield. 



D. Mechanisms of assembly 

Further information about the mechanisms of self- 
assembly can be gleaned from Fig. [5] which shows how 
the average cluster size evolves with MC time for dif- 
ferent temperatures at the optimal patch width. One 
striking feature is that for T < 0.14 ek^ 1 the average 
cluster size goes through a maximum before decreasing 
towards twelve. The pathway to cluster formation is 
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FIG. 5: (Colour Online) Some thermodynamic properties of 
a 12-particle system, (a) The heat capacity as a function of 
temperature for different values of a. Each line is labelled by 
the value of a (given in radians), (b) The free energy asso- 
ciated with the order parameter Q plotted at four different 
temperatures, all at a patch width of a — 0.45. Each line is 
labelled by the temperature, with T = 0.174 being the value 
of Tciust for this patch width, (c) Equilibrium probability of 
Q = 12 (solid lines). For comparison the yield of icosahedra 
obtained after 80 000 MC cycles in self assembly simulations 
(dotted lines) is included. Each line is labelled with the value 
of a (given in radians). 
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FIG. 6: (Colour Online) The mean cluster size as a function 
of number of MC cycles at a = 0.45 and at different tempera- 
tures (as labelled) for a system of 1200 particles at a number 
density of 0.15<TjJj. Each line is an average over 10 simula- 
tions. 

through larger disordered clusters, which then "bud off" 
icosahedra, rather than through the growth of icosahe- 
dra from smaller units. When an icosahedron with all 
the patches correctly oriented forms within a larger clus- 
ter, it will only interact very weakly with the rest of the 
cluster, and so will be able to escape relatively easily. 
The height of the maximum in the average cluster size 
increases with decreasing temperature because the rate 
of cluster growth increases compared to the rate at which 
icosahedra are annealed out. 

These results indicate therefore that the self-assembly 
of icosahedra is facilitated by the formation of a 
metastable intermediate phase, and the liquid phase, 
which we previously noted is stable at larger cr, plays an 
important role, even when it has disappeared from the 
equilibrium phase diagram. This system therefore pro- 
vides an example of Ostwald's step rule^ where nucle- 
ation initially leads to a metastable phase. For example, 
protein crystallization can be enhanced by the formation 
of a metastable protein-rich phase from which the crystal 
can more easily nucleate! 53 ' 54 

Fig.[7]illustrates some of the changes in assembly mech- 
anism as a function of temperature and patch width. In 
general, conditions of high temperature and low a favour 
a mechanism in which capsids are built up by a stepwise 
addition of monomers and small clusters. In Fig. [7Ja) 
the cluster size distribution is dominated by monomers 
and small intermediates after 80 MC cycles. Over the 
course of the simulation the population of particles shifts 
steadily into icosahedra, with almost no clusters larger 
than 12 particles being formed. Conversely, conditions 
of low temperature and high a (close to the range of pa- 
rameters where there is a stable liquid phase 5 - 5 -) favour the 
formation of large clusters, which then anneal to icosa- 
hedra. Fig. [TJb) shows the rapid formation of very large 
clusters at the start of the simulation. The average clus- 



ter size then shrinks back down by the "budding" mecha- 
nism mentioned above until the majority of particles are 
in icosahedra. It is interesting to note that at long times 
the large clusters that are left preferentially adopt sizes 
that are magic numbers for the isotropic Lennard- Jones 
potential, e.g. the peak at N = 19 in Fig.[7]Jb) that corre- 
sponds to two interpenetrating centred icosahedra4£ At 
the optimal conditions for assembly (Fig. [3c)), a combi- 
nation of the two mechanisms seems to operate. 

One effect of the availability of this budding mecha- 
nism is to increase the range of temperature over which 
icosahedra assemble. For example, it is noticeable from 
Fig. [3c) that as a decreases and the system moves fur- 
ther from the region of parameter space where the liquid 
phase is stable, the temperature range for which there 
are significant yields decreases significantly. 

The budding mechanism of self-assembly is very dif- 
ferent from that so far observed in experimenta l 21 ' 26 ' 56 
and theoretical (both kinetic model s 16 ' 56 and direct 
simulations 15 ) work on viruses, where direct nucleation 
of the virus capsid by stepwise addition of individual cap- 
somers is the norm. We suspect that the greater speci- 
ficity of the interactions in real viruses both makes the 
formation of large aggregates less likely under conditions 
where the complete capsids are most stable, and prevents 
the easy rearrangement of subunit particles in any large 
aggregates that form. 

E. Dependence on density 

All of the results described thus far are for simulations 
at the same density, 0.15er£"j, which raises the question 
of how the assembly behaviour depends on particle den- 
sity. Fig. [5fa) shows the yield of icosahedra at a range of 
particle densities, for the value of a that we found to be 
optimal for a density of 0.15 a£j 3 . The dependence of the 
yield on density is weak, with markedly decreased yields 
only observed at very low densities, where the times re- 
quired for the particles to diffuse and reach each other are 
significantly longer. At very low densities, assembly also 
becomes comparatively more successful at low tempera- 
tures. This small effect arises because the formation of 
large aggregates, which generally competes with the as- 
sembly of icosahedra at lower temperature, is disfavoured 
at low density (Fig.[5fb)), and so the successful growth of 
icosahedra by monomer addition is somewhat enhanced. 



F. Effect of inaccurate patch placement 

If self-assembling structures are to be constructed from 
nanoparticles and colloids, experimental methods will 
have to be developed to pattern the surfaces of the par- 
ticles appropriately with "sticky" and "non-sticky" re- 
gions. Since a certain degree of variability is intrinsic to 
synthetic colloids, it is useful to consider how robust the 
self-assembly process is with respect to imperfect place- 
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FIG. 7: (Colour Online) The fraction of particles in clusters 
of a certain size as a function of cluster size at different sim- 
ulation times for parameters (a) close to the transition to a 
monomer gas, T = 0.16, a — 0.4 (b) close to the transition to 
a liquid, T = 0.14, a = 0.55 and (c) giving optimal conditions 
for assembly, T — 0.14, a = 0.45. In each case the data is 
an average over 100 simulations, where each simulation is of 
1200 particles at a number density of 0.15 <t£j. 



merit of the patches. These factors may also be impor- 
tant in the consideration of the evolution of proteins that 
self-assemble into some quaternary structure, as the ro- 
bustness of self-assembly will have a strong effect on how 
easily such proteins can evolve. 



0.3 
0.28 
0.26 
0.24 
0.22 

0.2 
0.18 
0.16 
0.14 
0.12 

0.1 
0.08 
0.06 
0.04 
0.02 




0.05 0.1 



0.15 0.2 0.25 
Density / o LJ ~ 3 



0.3 0.35 0.4 



(b) 




0.15 0.2 0.25 
Density / Ojj" 



FIG. 8: (Colour Online) (a) The number of icosahedra formed 
and (b) the mean cluster size (averaged over particles) after 
80 000 MC cycles as a function of density and temperature for 
1200 particles at a fixed patch width of a = 0.45. 



In order to examine the effect of imperfect patch place- 
ment on assembly, simulations were performed where a 
random noise term was added to the patch vectors for 
each particle, with the noise for each patch on each par- 
ticle being generated separately. The modified patch vec- 
tors p are given by 



Po + Mgi 
|pb+MPi| 



(7) 



where p~o is the original (noiseless) patch vector, p! is 
a random unit vector, and [i is a parameter that deter- 
mines the magnitude of the noise in the patch placement. 
The average angle of deviation of p from p"o is approx- 
imately linear in /Lt, with a proportionality constant of 
0.78 radians. 

Analogous plots to Fig. G2a) are depicted in Fig. [9] for 
increasing values of the parameter (i. A low level of noise 
(fi = 0.1) seems to have little effect on the formation 
of icosahedra. Yields away from the optimum parame- 
ters are generally slightly depressed, and the temperature 
above which icosahedra are not observed is slightly re- 
duced. With increasing noise levels however, the region 
in the (<r, T) plane where assembly is possible shrinks 
rapidly, and by [i = 0.3 it has nearly disappeared. The 
noise has the general effect of reducing the binding energy 
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FIG. 9: (Colour Online) The effect of noisy patch placement. 
The heat maps show the numbers of icosahedra formed after 
80 000 MC cycles as a function of the patch width a (mea- 
sured in radians) and the temperature for 1200 particles at 
a number density of 0.15 <7jjj for different noise levels in the 
patch placement: (a) fi — 0.1, (b) /i = 0.2 and (c) fi = 0.3. 



of an icosahedral structure. Thus, the transition temper- 
ature to a monomer gas is reduced, as is the patch width 
above which larger aggregates become competitive. At 
low a, i.e. narrow patches, the reduction in the binding 
energy is more pronounced, and icosahedra again become 
unstable. 

The maximum yield of icosahedra as a function of fi 
is shown in Fig. 1101 Interestingly, the maximum yield 
is virtually independent of noise level up to fi = 0.2, 
after which it falls off rapidly. Beyond fj, — 0.35 there is 
virtually no successful assembly, /i = 0.2 corresponds to 
a mean deviation of 0.16 radians, or around 9°, and may 
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FIG. 10: The maximum yield of icosahedra as a function of 
the patch placement noise parameter /j,. In each case the 
yield is the average of five simulations at the values of T and 
a which give the maximum yield for that value of fi. Each 
simulation contains 1200 particles, at a number density of 
0.15 a" 3 . 



represent a target for the accuracy of patch placement in 
experiments, although one would expect the threshold to 
be somewhat dependent on the target structure. 



G. Use of a single ring patch 

The experimental production of five distinct patches 
that are accurately placed on a nanoparticle or colloidal 
surface is likely to be very challenging. One potentially 
easier alternative may be the creation of a continuous 
attractive ring on the particle's surfaced To mimic the 
effect of the discrete patches the ring should pass through 
the positions at which the individual patches would oth- 
erwise have been placed. 

To describe such a case the angular modulation of the 
interaction potential in the model (Eq. 2]) is replaced by 



Gij (iij,fli) = exp - 



2ct 2 



(8) 



where the orientation of the ring patch is described by a 
patch vector that passes through the centre of the ring, 
v is the cone angle from the patch vector at which the 
attraction is greatest, and a is a measure of the angular 
width of the ring. The yields of icosahedra for particles 
with a single ring patch at the appropriate value of v, 
namely v — 1.017 radians, are shown in Fig. 1111 

A number of comparisons can be made to the simula- 
tions where five separate patches were used. Firstly, the 
maximum yield observed when using a ring patch is the 
same to within statistical error (85% as opposed to 88% 
obtained with five discrete patches). It seems that the 
loss of specificity in changing to a ring patch is not too 
important, since the value of v is chosen such that the 
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FIG. 11: (Colour Online) Icosahedron yields using ring 
patches. The heat map shows the number of icosahedra 
formed after 80 000 MC cycles as a function of the patch width 
a (measured in radians) and the temperature for 1200 parti- 
cles at a number density of 0.15 a^j , using a single ring patch 
on each particle with v = 1.017 radians. 



most favourable configuration is still one where each of 
the particles has five neighbours. Secondly, the icosahe- 
dra are also stable to a slightly higher temperature when 
using ring patches, because each particle in an icosahe- 
dron is constrained in only two orientational degrees of 
freedom rather than three, reducing the entropy differ- 
ence between the icosahedra and the monomeric vapour. 
Thirdly, the optimal region is shifted to lower a, because 
for a given value of a more of the surface of the parti- 
cles is attractive than for particles with discrete patches 
at the same value of a. Consequently liquid-like aggre- 
gates become competitive at lower values of a. Finally, 
the assembly process is somewhat less vulnerable to ki- 
netic traps, because the continuous ring patch facilitates 
rearrangement of the particles. 

Our results suggest that ring patches represent an at- 
tractive alternative method of patterning particle sur- 
faces to form specific target structures, with prospects 
for easier synthesis and comparable yields to those ob- 
tained using discrete patches. 



H. Effect of a short-range potential 

All of the previous simulations have used an orienta- 
tionally modulated version of the Lennard-Jones poten- 
tial given in Eq. Q] However, the interaction potentials 
for real colloids and proteins tend to be much shorter 
in ranged Using a longer range potential has computa- 
tional advantages, but in order to investigate the effect of 
the range of the potential, a set of simulations was per- 
formed using a generalised Lennard-Jones potential with 
exponents of 20 and 10: 



f L T'»(r)=4 f 



(fH)-_(£H) 
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FIG. 12: (Colour Online) (a) The number of icosahedra 
formed and (b) the mean cluster size (averaged over particles) 
after 80 000 MC cycles using a shorter-range 20-10 potential 
as a function of the patch width a and temperature for 1200 
particles at a number density of 0.15 a^f. 



The results of these simulations are shown in Fig. [121 

The overall trends are similar to those for the longer 
range potential. However, yields are generally reduced, 
because the shorter range potential makes it less likely for 
particles to enter each others' attractive potential wells 
and form bonds, and because the barriers to rearrange- 
ments increase.— As a result the entire self assembly pro- 
cess proceeds more slowly. At longer simulation times the 
yields recover, with a maximum yield of icosahedron of 
97% after 1.2 x 10 6 MC cycles. Another general effect 
is that because of the slower bond formation, any ag- 
gregates that form will tend to be reduced in size (Fig. 

TO). 

The shorter range of the 20-10 potential energeti- 
cally destabilises the liquid phase ) 59 i 60 thus increasing 
the maximum value of a at which the system is stable as 
a collection of icosahedra. Hence, icosahedra are able to 
assemble up to higher a. Furthermore, since the budding 
mechanism of assembly requires the formation of a liq- 
uid droplet, assembly by nucleation is also comparatively 
more dominant than in the case of the 12-6 potential. 

Decreasing the range of the potential also has some 
other small effects, seen in Fig. lT27 a). The vibrational en- 
tropy of the icosahedra is reduced, lowering T c i us t slightly. 
There is also a small increase in the value of a below 
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which the kinetic suppression of icosahedron formation 
occurs, since at low a the short range further increases 
the difficulty of assembly. 

IV. CONCLUSIONS 

By making use of a minimal model, we have been able 
to simulate the behaviour of patchy particles designed 
to assemble into monodisperse icosahedral clusters. We 
found that, given the right conditions, the assembly pro- 
cess was rapid, efficient and robust. By performing ar- 
rays of simulations at different temperatures and using 
different patch widths, we were able to map out the re- 
gions of parameter space in which assembly is successful, 
and to identify the major competing behaviours arising 
in the regions of parameter space which did not lead to 
successful assembly. 

From these results we can learn something about the 
general principles required when designing objects to self- 
assemble. The optimal patch width represents a compro- 
mise between the energetic and kinetic requirements for 
self-assembly. As the patch width is decreased, the en- 
ergy gap between the target structure and other possible 
competing structures increases. By contrast, kinetic ac- 
cessibility improves as a increases, because the system 
is able to escape from incorrect configurations more eas- 
ily. The lesson is that the interactions need to be specific 
enough to sufficiently favour the target structure, but 
that overspecifying them can inhibit the dynamics of as- 
sembly. 

We also found that the dynamic pathways to self- 
assembly can be complex and non-intuitive. We ob- 
served, for example, that, besides the better known nu- 
cleation pathways, clusters could also form through disor- 
dered intermediates. These mechanisms are less likely to 
be important in biological systems, but may be relevant 
for synthetic self-assembling systems. In our simulations 
the optimal yields were found under conditions where 
both mechanisms made contributions to the assembly. 

Assembly via nucleation pathways appears to be much 
more successful close to the maximum temperature at 
which icosahedra are stable. Kinetic traps are avoided, 
and alternative processes such as the formation of liq- 
uid droplets are inhibited. Since high temperatures 
are equivalent to weak bonding interactions, this find- 
ing is consistent with the experimental observation that 
protein-protein interactions in virus capsids are generally 
weakj^i 



The assembly process is not strongly affected by the 
particle density. Changing the range of the potential 
also has little effect on the yields at long times (although 
assembly mechanisms involving liquid droplets are inhib- 
ited by a shorter range potential). The relative indepen- 
dence of the behaviour of the system with respect to these 
factors indicates that our findings are quite general, and 
not just specific to the parameter range we focussed on. 

The self-assembly process is encouragingly robust with 
respect to the design and placement of patches. Simula- 
tions in which the patches were positioned with random 
inaccuracies still produced almost unchanged yields of 
correctly formed icosahedra up to a threshold level of 
average error. While most of our simulations used par- 
ticles with five discrete patches, we found that particles 
with a single ring-shaped patch, which may be easier to 
produce experimentally, assembled with a similar level of 
efficiency. Both of these findings may be seen as very en- 
couraging for experimentalists seeking to synthesise self- 
assembling nanoparticles and colloids. 

Our study focussed on a simple model in order to be 
able to explore the design space in detail. Although the 
patchy spheres of our model clearly lack the complexity 
of the units involved in monodisperse self-assembly in bi- 
ological systems, the comprehensive understanding that 
this simplicity allows us to achieve should provide useful 
insights into biological examples. Furthermore, this fea- 
ture may be an advantage when considering the design 
of synthetic self-assembling building blocks. 

The target structure used in this work was relatively 
simple, and formation of monodisperse clusters was cor- 
respondingly rapid. In future work we intend to explore 
how the ease of assembly depends on the geometry of the 
target structure, and to investigate the limits on the com- 
plexity of structures that can be generated from this rel- 
atively simple interaction scheme. In addition, there are 
various features which could be added to the potential, 
most notably a dependence on torsional angles. We ex- 
pect that a torsional dependence would greatly increase 
the range of achievable target structures, and also pro- 
duce behaviour more closely mirroring that of biological 
systems. 
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